#############################################################
### Replication Code: Main Analysis                       ###
### Title: Can Conservatives Be Persuaded to Support UBI? ###
### Author: Eddy S. F. Yeung                              ###
### Version: September 6, 2022                            ###
#############################################################

# Set-up ----
## Clean the R environment and set the working directory
rm(list = ls())
setwd("~/Desktop/UBI/PB_replication") # change to your own directory

## Load the required packages
library(tidyverse) # version 1.3.1
library(multcomp)  # version 1.4-17
library(broom)     # version 1.0.0
library(estimatr)  # version 0.30.4
library(texreg)    # version 1.37.5
library(gridExtra) # version 2.3
library(extrafont) # version 0.17

## Import the cleaned dataset
load("cleaned.RData")

# Analysis ----
## Summary statistics for conservatives' UBI support across experimental groups ----
### Mean
with(df, c(mean(support[group == 0 & conservative == 1]), 
           mean(support[group == 1 & conservative == 1]),
           mean(support[group == 2 & conservative == 1])))

### Sample size
with(df, c(length(support[group == 0 & conservative == 1]), 
           length(support[group == 1 & conservative == 1]),
           length(support[group == 2 & conservative == 1])))

### Standard deviation
with(df, c(sd(support[group == 0 & conservative == 1]), 
           sd(support[group == 1 & conservative == 1]),
           sd(support[group == 2 & conservative == 1])))

## Table 1: OLS Regression Corresponding to Equation 1 ----
### Specification 1
mod1 <- lm_robust(support ~ treatment1 + treatment2 + conservative +
                    treatment1 * conservative + treatment2 * conservative,
                  data = df)

### Specification 2
mod2 <- lm_robust(support ~ treatment1 + treatment2 + conservative +
                    treatment1 * conservative + treatment2 * conservative +
                    age + female + married + educ,
                  data = df)

### Specification 3
mod3 <- lm_robust(support ~ treatment1 + treatment2 + conservative +
                    treatment1 * conservative + treatment2 * conservative +
                    age + female + married + educ + income + white + black +
                    employ + union + pid + welfare.receive,
                  data = df)

### Specification 4
mod4 <- lm_robust(support ~ treatment1 + treatment2 + conservative +
                    treatment1 * conservative + treatment2 * conservative +
                    age + female + married + educ + income + white + black +
                    employ + union + pid + welfare.receive + welfare.bad +
                    racial.resent + demo.change,
                  data = df)

### Export the regression table to LaTeX
texreg(list(mod1, mod2, mod3, mod4),
       stars = c(0.01, 0.05, 0.10),
       include.ci = F,
       custom.header = list("Dependent Variable: Seven-Point UBI Support" = 1:4),
       custom.coef.names = 
         c("Constant", "Equalizing-Opportunity (EO) Frame",
           "Limiting-Government (LG) Frame", "Conservative (= 1)",
           "EO Frame × Conservative", "LG Frame × Conservative",
           "Age", "Female (= 1)", "Married (= 1)", "Education (7 = Highest)", 
           "Household Income (12 = Highest)", "White (= 1)", "Black (= 1)", 
           "Not Employed (= 1)", "Retired (= 1)", "Household Union Membership (= 1)", 
           "Party ID (6 = Strong Republican)", "Welfare Recipient (= 1)", 
           "Welfare System Evaluation (1 = Bad)", "Racial Resentment (8 = Highest)", 
           "Demographic Change (1 = Yes)"),
       custom.note = "Entries are OLS estimates with robust standard errors in parentheses.
       All significance tests are two-tailed with the following notations:
       $^{***}p<0.01$; $^{**}p<0.05$; $^{*}p<0.1$.",
       caption = "OLS Regression Corresponding to Equation 1",
       fontsize = "small")

## H1b: test whether beta4 is greater than zero ----
temp.t <- coef(summary(mod4))["treatment1:conservative", "t value"] # extract t-value
temp.df <- mod4$df["treatment1:conservative"] # extract degrees of freedom
pt(q = temp.t, df = temp.df, lower.tail = F) # one-sided p-value

## H2b: test whether beta5 is greater than zero ----
temp.t <- coef(summary(mod4))["treatment2:conservative", "t value"] # extract t-value
temp.df <- mod4$df # extract degrees of freedom
pt(q = temp.t, df = temp.df, lower.tail = F) # one-sided p-value

## Figure 2: Average Support for UBI among Conservatives across Experimental Groups ----
summary_df <-
  df[which(df$conservative == 1), ] %>%
  group_by(group) %>%
  do(tidy(lm_robust(support ~ 1, data = .))) %>%
  mutate(support = estimate)

fig_mean.support <-
  ggplot(summary_df, aes(x = group, y = support)) +
  geom_point(data = df[which(df$conservative == 1), ],
             position = position_jitter(width = 0.1, height = 0.2, seed = 1234567),
             alpha = 0.2) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), size = .9, width = 0) +
  theme_bw() +
  scale_x_discrete(labels = c("0" = "Control Group\n(No Frame)",
                              "1" = "Treatment Group 1\n(Equalizing Opportunity)",
                              "2" = "Treatment Group 2\n(Limiting Government)")) +
  scale_y_continuous(breaks = seq(0, 6, by = 1)) +
  xlab("") +
  ylab("UBI Support among Conservatives\n(0 = Strongly Against; 6 = Strongly In Favor)") +
  coord_cartesian(ylim = c(-0.1, 6.1)) +
  theme(text = element_text(family = "Times", size = 14),
        axis.text = element_text(family = "Times", size = 13))
ggsave(file = "fig_mean.support.pdf", fig_mean.support, width = 7.5, height = 5)

## Figure 3: Conditional Average Treatment Effects on Conservatives ----
### Create an empty data frame to store the estimates for later visualization
df.CATE1 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE1) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE1[1:3, 1] <- 1
df.CATE1[4:6, 1] <- 2

### H1a: Test whether (beta1 + beta4) is different from zero (7-point scale DV)
### All respondents (including inattentive respondents)
H1a.full <- glht(mod3, linfct = c("treatment1 + treatment1:conservative = 0"))
df.CATE1[1, 2] <- confint(H1a.full)$confint[ , "Estimate"]
df.CATE1[1, 3] <- confint(H1a.full)$confint[ , "lwr"]
df.CATE1[1, 4] <- confint(H1a.full)$confint[ , "upr"]
df.CATE1[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3.alt1 <- lm_robust(support ~ treatment1 + treatment2 + conservative +
                         treatment1 * conservative + treatment2 * conservative +
                         age + female + married + educ + income + white + black +
                         employ + union + pid + welfare.receive,
                       data = subset(df, attention == 1 | attention == 2))
H1a.alt1 <- glht(mod3.alt1, linfct = c("treatment1 + treatment1:conservative = 0"))
df.CATE1[2, 2] <- confint(H1a.alt1)$confint[ , "Estimate"]
df.CATE1[2, 3] <- confint(H1a.alt1)$confint[ , "lwr"]
df.CATE1[2, 4] <- confint(H1a.alt1)$confint[ , "upr"]
df.CATE1[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3.alt2 <- lm_robust(support ~ treatment1 + treatment2 + conservative +
                         treatment1 * conservative + treatment2 * conservative +
                         age + female + married + educ + income + white + black +
                         employ + union + pid + welfare.receive,
                       data = subset(df, attention == 2))
H1a.alt2 <- glht(mod3.alt2, linfct = c("treatment1 + treatment1:conservative = 0"))
df.CATE1[3, 2] <- confint(H1a.alt2)$confint[ , "Estimate"]
df.CATE1[3, 3] <- confint(H1a.alt2)$confint[ , "lwr"]
df.CATE1[3, 4] <- confint(H1a.alt2)$confint[ , "upr"]
df.CATE1[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (7-point scale DV)
### All respondents (including inattentive respondents)
H2a.full <- glht(mod3, linfct = c("treatment2 + treatment2:conservative = 0"))
df.CATE1[4, 2] <- confint(H2a.full)$confint[ , "Estimate"]
df.CATE1[4, 3] <- confint(H2a.full)$confint[ , "lwr"]
df.CATE1[4, 4] <- confint(H2a.full)$confint[ , "upr"]
df.CATE1[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.alt1 <- glht(mod3.alt1, linfct = c("treatment2 + treatment2:conservative = 0"))
df.CATE1[5, 2] <- confint(H2a.alt1)$confint[ , "Estimate"]
df.CATE1[5, 3] <- confint(H2a.alt1)$confint[ , "lwr"]
df.CATE1[5, 4] <- confint(H2a.alt1)$confint[ , "upr"]
df.CATE1[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.alt2 <- glht(mod3.alt2, linfct = c("treatment2 + treatment2:conservative = 0"))
df.CATE1[6, 2] <- confint(H2a.alt2)$confint[ , "Estimate"]
df.CATE1[6, 3] <- confint(H2a.alt2)$confint[ , "lwr"]
df.CATE1[6, 4] <- confint(H2a.alt2)$confint[ , "upr"]
df.CATE1[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE1$group <- factor(df.CATE1$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                    "Treatment Group 2\n(Limiting Government)"))
fig_CATE1 <- 
  ggplot(data = df.CATE1, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr), position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Conservatives") +
  ggtitle("Panel A: Seven-Point UBI Support") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.title = element_blank(),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-1.15, 1.15)) +
  guides(color = guide_legend(nrow = 1, byrow = T))
fig_CATE1

### Create an empty data frame to store the estimates for later visualization
df.CATE2 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE2) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE2[1:3, 1] <- 1
df.CATE2[4:6, 1] <- 2

### H1a: Test whether (beta1 + beta4) is different from zero (binary DV)
### All respondents (including inattentive respondents)
mod3.binary0 <- lm_robust(percent.sup ~ treatment1 + treatment2 + conservative +
                            treatment1 * conservative + treatment2 * conservative +
                            age + female + married + educ + income + white + black +
                            employ + union + pid + welfare.receive,
                          data = df)
H1a.binary0 <- glht(mod3.binary0, linfct = c("treatment1 + treatment1:conservative = 0"))
df.CATE2[1, 2] <- confint(H1a.binary0)$confint[ , "Estimate"]
df.CATE2[1, 3] <- confint(H1a.binary0)$confint[ , "lwr"]
df.CATE2[1, 4] <- confint(H1a.binary0)$confint[ , "upr"]
df.CATE2[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3.binary1 <- lm_robust(percent.sup ~ treatment1 + treatment2 + conservative +
                            treatment1 * conservative + treatment2 * conservative +
                            age + female + married + educ + income + white + black +
                            employ + union + pid + welfare.receive,
                          data = subset(df, attention == 1 | attention == 2))
H1a.binary1 <- glht(mod3.binary1, linfct = c("treatment1 + treatment1:conservative = 0"))
df.CATE2[2, 2] <- confint(H1a.binary1)$confint[ , "Estimate"]
df.CATE2[2, 3] <- confint(H1a.binary1)$confint[ , "lwr"]
df.CATE2[2, 4] <- confint(H1a.binary1)$confint[ , "upr"]
df.CATE2[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3.binary2 <- lm_robust(percent.sup ~ treatment1 + treatment2 + conservative +
                            treatment1 * conservative + treatment2 * conservative +
                            age + female + married + educ + income + white + black +
                            employ + union + pid + welfare.receive,
                          data = subset(df, attention == 2))
H1a.binary2 <- glht(mod3.binary2, linfct = c("treatment1 + treatment1:conservative = 0"))
df.CATE2[3, 2] <- confint(H1a.binary2)$confint[ , "Estimate"]
df.CATE2[3, 3] <- confint(H1a.binary2)$confint[ , "lwr"]
df.CATE2[3, 4] <- confint(H1a.binary2)$confint[ , "upr"]
df.CATE2[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (binary DV)
### All respondents (including inattentive respondents)
H2a.binary0 <- glht(mod3.binary0, linfct = c("treatment2 + treatment2:conservative = 0"))
df.CATE2[4, 2] <- confint(H2a.binary0)$confint[ , "Estimate"]
df.CATE2[4, 3] <- confint(H2a.binary0)$confint[ , "lwr"]
df.CATE2[4, 4] <- confint(H2a.binary0)$confint[ , "upr"]
df.CATE2[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.binary1 <- glht(mod3.binary1, linfct = c("treatment2 + treatment2:conservative = 0"))
df.CATE2[5, 2] <- confint(H2a.binary1)$confint[ , "Estimate"]
df.CATE2[5, 3] <- confint(H2a.binary1)$confint[ , "lwr"]
df.CATE2[5, 4] <- confint(H2a.binary1)$confint[ , "upr"]
df.CATE2[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.binary2 <- glht(mod3.binary2, linfct = c("treatment2 + treatment2:conservative = 0"))
df.CATE2[6, 2] <- confint(H2a.binary2)$confint[ , "Estimate"]
df.CATE2[6, 3] <- confint(H2a.binary2)$confint[ , "lwr"]
df.CATE2[6, 4] <- confint(H2a.binary2)$confint[ , "upr"]
df.CATE2[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE2$group <- factor(df.CATE2$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                   "Treatment Group 2\n(Limiting Government)"))
fig_CATE2 <- 
  ggplot(data = df.CATE2, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr), position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Conservatives") +
  ggtitle("Panel B: Percentage Support for UBI") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-.3, .3)) +
  guides(colour = guide_legend(nrow = 1, byrow = TRUE))
fig_CATE2

### Combine the two plots (fig_CATE1 and fig_CATE2) and let them share a common legend
get_legend <- function(myggplot) {
  temp <- ggplot_gtable(ggplot_build(myggplot))
  leg <- which(sapply(temp$grobs, function(x) x$name) == "guide-box")
  legend <- temp$grobs[[leg]]
  return(legend)
}
legend <- get_legend(fig_CATE1) # save the legend
fig_CATE <- 
  grid.arrange(arrangeGrob(fig_CATE1 + theme(legend.position = "none"),
                           fig_CATE2 + theme(legend.position = "none"),
                           nrow = 1),
               legend, nrow = 2, heights = c(10, 1))
ggsave(file = "fig_CATE.pdf", fig_CATE, width = 8, height = 4.5)

## Figure 4: Conditional Average Treatment Effects on Conservatives after Redefining the Conservatism Score Threshold ----
### Create an empty data frame to store the estimates for later visualization
df.CATE3 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE3) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE3[1:3, 1] <- 1
df.CATE3[4:6, 1] <- 2

### H1a: Test whether (beta1 + beta4) is different from zero (cons.score >= 6) 
### All respondents (including inattentive respondents)
mod3.con6.1 <- lm_robust(support ~ treatment1 + treatment2 + conservative6 +
                           treatment1 * conservative6 + treatment2 * conservative6 +
                           age + female + married + educ + income + white + black +
                           employ + union + pid + welfare.receive,
                         data = df)
H1a.con6.1 <- glht(mod3.con6.1, linfct = c("treatment1 + treatment1:conservative6 = 0"))
df.CATE3[1, 2] <- confint(H1a.con6.1)$confint[ , "Estimate"]
df.CATE3[1, 3] <- confint(H1a.con6.1)$confint[ , "lwr"]
df.CATE3[1, 4] <- confint(H1a.con6.1)$confint[ , "upr"]
df.CATE3[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3.con6.2 <- lm_robust(support ~ treatment1 + treatment2 + conservative6 +
                           treatment1 * conservative6 + treatment2 * conservative6 +
                           age + female + married + educ + income + white + black +
                           employ + union + pid + welfare.receive,
                         data = subset(df, attention == 1 | attention == 2))
H1a.con6.2 <- glht(mod3.con6.2, linfct = c("treatment1 + treatment1:conservative6 = 0"))
df.CATE3[2, 2] <- confint(H1a.con6.2)$confint[ , "Estimate"]
df.CATE3[2, 3] <- confint(H1a.con6.2)$confint[ , "lwr"]
df.CATE3[2, 4] <- confint(H1a.con6.2)$confint[ , "upr"]
df.CATE3[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3.con6.3 <- lm_robust(support ~ treatment1 + treatment2 + conservative6 +
                           treatment1 * conservative6 + treatment2 * conservative6 +
                           age + female + married + educ + income + white + black +
                           employ + union + pid + welfare.receive,
                         data = subset(df, attention == 2))
H1a.con6.3 <- glht(mod3.con6.3, linfct = c("treatment1 + treatment1:conservative6 = 0"))
df.CATE3[3, 2] <- confint(H1a.con6.3)$confint[ , "Estimate"]
df.CATE3[3, 3] <- confint(H1a.con6.3)$confint[ , "lwr"]
df.CATE3[3, 4] <- confint(H1a.con6.3)$confint[ , "upr"]
df.CATE3[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (cons.score >= 6) 
### All respondents (including inattentive respondents)
H2a.con6.1 <- glht(mod3.con6.1, linfct = c("treatment2 + treatment2:conservative6 = 0"))
df.CATE3[4, 2] <- confint(H2a.con6.1)$confint[ , "Estimate"]
df.CATE3[4, 3] <- confint(H2a.con6.1)$confint[ , "lwr"]
df.CATE3[4, 4] <- confint(H2a.con6.1)$confint[ , "upr"]
df.CATE3[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.con6.2 <- glht(mod3.con6.2, linfct = c("treatment2 + treatment2:conservative6 = 0"))
df.CATE3[5, 2] <- confint(H2a.con6.2)$confint[ , "Estimate"]
df.CATE3[5, 3] <- confint(H2a.con6.2)$confint[ , "lwr"]
df.CATE3[5, 4] <- confint(H2a.con6.2)$confint[ , "upr"]
df.CATE3[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.con6.3 <- glht(mod3.con6.3, linfct = c("treatment2 + treatment2:conservative6 = 0"))
df.CATE3[6, 2] <- confint(H2a.con6.3)$confint[ , "Estimate"]
df.CATE3[6, 3] <- confint(H2a.con6.3)$confint[ , "lwr"]
df.CATE3[6, 4] <- confint(H2a.con6.3)$confint[ , "upr"]
df.CATE3[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE3$group <- factor(df.CATE3$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                    "Treatment Group 2\n(Limiting Government)"))
fig_CATE3 <- 
  ggplot(data = df.CATE3, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr), position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Conservatives") +
  ggtitle("Panel A: Conservatism Score Threshold = 6") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.title = element_blank(),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-1.15, 1.15)) +
  guides(color = guide_legend(nrow = 1, byrow = T))
fig_CATE3

### Create an empty data frame to store the estimates for later visualization
df.CATE4 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE4) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE4[1:3, 1] <- 1
df.CATE4[4:6, 1] <- 2

### H2a: Test whether (beta1 + beta4) is different from zero (cons.score >= 7)
### All respondents (including inattentive respondents)
mod3.con7.1 <- lm_robust(support ~ treatment1 + treatment2 + conservative7 +
                           treatment1 * conservative7 + treatment2 * conservative7 +
                           age + female + married + educ + income + white + black +
                           employ + union + pid + welfare.receive,
                         data = df)
H1a.con7.1 <- glht(mod3.con7.1, linfct = c("treatment1 + treatment1:conservative7 = 0"))
df.CATE4[1, 2] <- confint(H1a.con7.1)$confint[ , "Estimate"]
df.CATE4[1, 3] <- confint(H1a.con7.1)$confint[ , "lwr"]
df.CATE4[1, 4] <- confint(H1a.con7.1)$confint[ , "upr"]
df.CATE4[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3.con7.2 <- lm_robust(support ~ treatment1 + treatment2 + conservative7 +
                           treatment1 * conservative7 + treatment2 * conservative7 +
                           age + female + married + educ + income + white + black +
                           employ + union + pid + welfare.receive,
                         data = subset(df, attention == 1 | attention == 2))
H1a.con7.2 <- glht(mod3.con7.2, linfct = c("treatment1 + treatment1:conservative7 = 0"))
df.CATE4[2, 2] <- confint(H1a.con7.2)$confint[ , "Estimate"]
df.CATE4[2, 3] <- confint(H1a.con7.2)$confint[ , "lwr"]
df.CATE4[2, 4] <- confint(H1a.con7.2)$confint[ , "upr"]
df.CATE4[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3.con7.3 <- lm_robust(support ~ treatment1 + treatment2 + conservative7 +
                           treatment1 * conservative7 + treatment2 * conservative7 +
                           age + female + married + educ + income + white + black +
                           employ + union + pid + welfare.receive,
                         data = subset(df, attention == 2))
H1a.con7.3 <- glht(mod3.con7.3, linfct = c("treatment1 + treatment1:conservative7 = 0"))
df.CATE4[3, 2] <- confint(H1a.con7.3)$confint[ , "Estimate"]
df.CATE4[3, 3] <- confint(H1a.con7.3)$confint[ , "lwr"]
df.CATE4[3, 4] <- confint(H1a.con7.3)$confint[ , "upr"]
df.CATE4[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (cons.score >= 7)
### All respondents (including inattentive respondents)
H2a.con7.1 <- glht(mod3.con7.1, linfct = c("treatment2 + treatment2:conservative7 = 0"))
df.CATE4[4, 2] <- confint(H2a.con7.1)$confint[ , "Estimate"]
df.CATE4[4, 3] <- confint(H2a.con7.1)$confint[ , "lwr"]
df.CATE4[4, 4] <- confint(H2a.con7.1)$confint[ , "upr"]
df.CATE4[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.con7.2 <- glht(mod3.con7.2, linfct = c("treatment2 + treatment2:conservative7 = 0"))
df.CATE4[5, 2] <- confint(H2a.con7.2)$confint[ , "Estimate"]
df.CATE4[5, 3] <- confint(H2a.con7.2)$confint[ , "lwr"]
df.CATE4[5, 4] <- confint(H2a.con7.2)$confint[ , "upr"]
df.CATE4[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.con7.3 <- glht(mod3.con7.3, linfct = c("treatment2 + treatment2:conservative7 = 0"))
df.CATE4[6, 2] <- confint(H2a.con7.3)$confint[ , "Estimate"]
df.CATE4[6, 3] <- confint(H2a.con7.3)$confint[ , "lwr"]
df.CATE4[6, 4] <- confint(H2a.con7.3)$confint[ , "upr"]
df.CATE4[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE4$group <- factor(df.CATE4$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                    "Treatment Group 2\n(Limiting Government)"))
fig_CATE4 <- 
  ggplot(data = df.CATE4, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr), position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Conservatives") +
  ggtitle("Panel B: Conservatism Score Threshold = 7") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "none",
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-1.15, 1.15)) +
  guides(colour = guide_legend(nrow = 1, byrow = TRUE))
fig_CATE4

### Combine the two plots (fig_CATE3 and fig_CATE4) and let them share a common legend
legend <- get_legend(fig_CATE3) # save the legend
fig_CATE_alt <- 
  grid.arrange(arrangeGrob(fig_CATE3 + theme(legend.position = "none"),
                           fig_CATE4 + theme(legend.position = "none"),
                           nrow = 1),
               legend, nrow = 2, heights = c(10, 1))
ggsave(file = "fig_CATE_alt.pdf", fig_CATE_alt, width = 8, height = 4.5)

## Figure 5: Conditional Average Treatment Effects on Economic and Social Conservatives ----
### Create an empty data frame to store the estimates for later visualization
df.CATE3 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE3) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE3[1:3, 1] <- 1
df.CATE3[4:6, 1] <- 2

### H1a: Test whether (beta1 + beta4) is different from zero (e.conservative)
### All respondents (including inattentive respondents)
mod3_econ0 <- lm_robust(support ~ treatment1 + treatment2 + e.conservative +
                          treatment1 * e.conservative + treatment2 * e.conservative +
                          age + female + married + educ + income + white + black +
                          employ + union + pid + welfare.receive,
                        data = df)
H1a.full <- glht(mod3_econ0, linfct = c("treatment1 + treatment1:e.conservative = 0"))
df.CATE3[1, 2] <- confint(H1a.full)$confint[ , "Estimate"]
df.CATE3[1, 3] <- confint(H1a.full)$confint[ , "lwr"]
df.CATE3[1, 4] <- confint(H1a.full)$confint[ , "upr"]
df.CATE3[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3_econ1 <- lm_robust(support ~ treatment1 + treatment2 + e.conservative +
                          treatment1 * e.conservative + treatment2 * e.conservative +
                          age + female + married + educ + income + white + black +
                          employ + union + pid + welfare.receive,
                        data = subset(df, attention == 1 | attention == 2))
H1a.alt1 <- glht(mod3_econ1, linfct = c("treatment1 + treatment1:e.conservative = 0"))
df.CATE3[2, 2] <- confint(H1a.alt1)$confint[ , "Estimate"]
df.CATE3[2, 3] <- confint(H1a.alt1)$confint[ , "lwr"]
df.CATE3[2, 4] <- confint(H1a.alt1)$confint[ , "upr"]
df.CATE3[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3_econ2 <- lm_robust(support ~ treatment1 + treatment2 + e.conservative +
                          treatment1 * e.conservative + treatment2 * e.conservative +
                          age + female + married + educ + income + white + black +
                          employ + union + pid + welfare.receive,
                        data = subset(df, attention == 2))
H1a.alt2 <- glht(mod3_econ2, linfct = c("treatment1 + treatment1:e.conservative = 0"))
df.CATE3[3, 2] <- confint(H1a.alt2)$confint[ , "Estimate"]
df.CATE3[3, 3] <- confint(H1a.alt2)$confint[ , "lwr"]
df.CATE3[3, 4] <- confint(H1a.alt2)$confint[ , "upr"]
df.CATE3[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (e.conservative)
### All respondents (including unattentive respondents)
H2a.full <- glht(mod3_econ0, linfct = c("treatment2 + treatment2:e.conservative = 0"))
df.CATE3[4, 2] <- confint(H2a.full)$confint[ , "Estimate"]
df.CATE3[4, 3] <- confint(H2a.full)$confint[ , "lwr"]
df.CATE3[4, 4] <- confint(H2a.full)$confint[ , "upr"]
df.CATE3[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.alt1 <- glht(mod3_econ1, linfct = c("treatment2 + treatment2:e.conservative = 0"))
df.CATE3[5, 2] <- confint(H2a.alt1)$confint[ , "Estimate"]
df.CATE3[5, 3] <- confint(H2a.alt1)$confint[ , "lwr"]
df.CATE3[5, 4] <- confint(H2a.alt1)$confint[ , "upr"]
df.CATE3[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.alt2 <- glht(mod3_econ2, linfct = c("treatment2 + treatment2:e.conservative = 0"))
df.CATE3[6, 2] <- confint(H2a.alt2)$confint[ , "Estimate"]
df.CATE3[6, 3] <- confint(H2a.alt2)$confint[ , "lwr"]
df.CATE3[6, 4] <- confint(H2a.alt2)$confint[ , "upr"]
df.CATE3[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE3$group <- factor(df.CATE3$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                    "Treatment Group 2\n(Limiting Government)"))
fig_CATE3 <- 
  ggplot(data = df.CATE3, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr), position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Economic Conservatives") +
  ggtitle("Panel A: Economic Conservatives") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.title = element_blank(),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-1.15, 1.15)) +
  guides(color = guide_legend(nrow = 1, byrow = T))
fig_CATE3

### Create an empty data frame to store the estimates for later visualization
df.CATE4 <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.CATE4) <- c("group", "estimate", "lwr", "upr", "sample")
df.CATE4[1:3, 1] <- 1
df.CATE4[4:6, 1] <- 2

### H1a: Test whether (beta1 + beta4) is different from zero (s.conservative)
### All respondents (including inattentive respondents)
mod3_soc0 <- lm_robust(support ~ treatment1 + treatment2 + s.conservative +
                         treatment1 * s.conservative + treatment2 * s.conservative +
                         age + female + married + educ + income + white + black +
                         employ + union + pid + welfare.receive,
                        data = df)
H1a.full <- glht(mod3_soc0, linfct = c("treatment1 + treatment1:s.conservative = 0"))
df.CATE4[1, 2] <- confint(H1a.full)$confint[ , "Estimate"]
df.CATE4[1, 3] <- confint(H1a.full)$confint[ , "lwr"]
df.CATE4[1, 4] <- confint(H1a.full)$confint[ , "upr"]
df.CATE4[1, 5] <- "All Respondents"

### Fully and partly attentive respondents only
mod3_soc1 <- lm_robust(support ~ treatment1 + treatment2 + s.conservative +
                         treatment1 * s.conservative + treatment2 * s.conservative +
                         age + female + married + educ + income + white + black +
                         employ + union + pid + welfare.receive,
                       data = subset(df, attention == 1 | attention == 2))
H1a.alt1 <- glht(mod3_soc1, linfct = c("treatment1 + treatment1:s.conservative = 0"))
df.CATE4[2, 2] <- confint(H1a.alt1)$confint[ , "Estimate"]
df.CATE4[2, 3] <- confint(H1a.alt1)$confint[ , "lwr"]
df.CATE4[2, 4] <- confint(H1a.alt1)$confint[ , "upr"]
df.CATE4[2, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
mod3_soc2 <- lm_robust(support ~ treatment1 + treatment2 + s.conservative +
                         treatment1 * s.conservative + treatment2 * s.conservative +
                         age + female + married + educ + income + white + black +
                         employ + union + pid + welfare.receive,
                        data = subset(df, attention == 2))
H1a.alt2 <- glht(mod3_soc2, linfct = c("treatment1 + treatment1:s.conservative = 0"))
df.CATE4[3, 2] <- confint(H1a.alt2)$confint[ , "Estimate"]
df.CATE4[3, 3] <- confint(H1a.alt2)$confint[ , "lwr"]
df.CATE4[3, 4] <- confint(H1a.alt2)$confint[ , "upr"]
df.CATE4[3, 5] <- "Fully Attentive Respondents Only"

### H2a: Test whether (beta2 + beta5) is different from zero (s.conservative)
### All respondents (including inattentive respondents)
H2a.full <- glht(mod3_soc0, linfct = c("treatment2 + treatment2:s.conservative = 0"))
df.CATE4[4, 2] <- confint(H2a.full)$confint[ , "Estimate"]
df.CATE4[4, 3] <- confint(H2a.full)$confint[ , "lwr"]
df.CATE4[4, 4] <- confint(H2a.full)$confint[ , "upr"]
df.CATE4[4, 5] <- "All Respondents"

### Fully and partly attentive respondents only
H2a.alt1 <- glht(mod3_soc1, linfct = c("treatment2 + treatment2:s.conservative = 0"))
df.CATE4[5, 2] <- confint(H2a.alt1)$confint[ , "Estimate"]
df.CATE4[5, 3] <- confint(H2a.alt1)$confint[ , "lwr"]
df.CATE4[5, 4] <- confint(H2a.alt1)$confint[ , "upr"]
df.CATE4[5, 5] <- "Fully and Partly Attentive Respondents Only"

### Fully attentive respondents only
H2a.alt2 <- glht(mod3_soc2, linfct = c("treatment2 + treatment2:s.conservative = 0"))
df.CATE4[6, 2] <- confint(H2a.alt2)$confint[ , "Estimate"]
df.CATE4[6, 3] <- confint(H2a.alt2)$confint[ , "lwr"]
df.CATE4[6, 4] <- confint(H2a.alt2)$confint[ , "upr"]
df.CATE4[6, 5] <- "Fully Attentive Respondents Only"

### Visualize the results
df.CATE4$group <- factor(df.CATE4$group,
                         levels = c(1:2),
                         labels = c("Treatment Group 1\n(Equalizing Opportunity)",
                                    "Treatment Group 2\n(Limiting Government)"))
fig_CATE4 <- 
  ggplot(data = df.CATE4, aes(x = group, y = estimate, color = sample, shape = sample)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey30", "grey60")) +
  scale_shape_manual(values = c(19, 15, 17)) +
  geom_errorbar(width = 0, aes(ymin = lwr, ymax = upr), position = position_dodge(.5)) +
  xlab("") + 
  ylab("CATE on Social Conservatives") +
  ggtitle("Panel B: Social Conservatives") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 13, family = "Times"),
        legend.position = "bottom",
        legend.direction = "horizontal",
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.title = element_blank(),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm"),
        plot.title = element_text(hjust = 0.5, size = 14, family = "Times")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  coord_cartesian(ylim = c(-1.15, 1.15)) +
  guides(color = guide_legend(nrow = 1, byrow = T))
fig_CATE4

### Combine the two plots (fig_CATE3 and fig_CATE4) and let them share a common legend
legend <- get_legend(fig_CATE3) # save the legend
fig_CATE_2d <- 
  grid.arrange(arrangeGrob(fig_CATE3 + theme(legend.position = "none"),
                           fig_CATE4 + theme(legend.position = "none"),
                           nrow = 1),
               legend, nrow = 2, heights = c(10, 1))
ggsave(file = "fig_CATE_2D.pdf", fig_CATE_2d, width = 8, height = 4.5)

## Figure 7: Individual Ideology and Association of UBI and TANF with Politicians from the Democratic Party ----
### Create an empty data frame to store the estimates for later visualization
df.asso <- data.frame(matrix(NA, nrow = 6, ncol = 5))
colnames(df.asso) <- c("sample", "estimate", "lwr", "upr", "Policy")
df.asso[1:2, 1] <- 1
df.asso[3:4, 1] <- 2
df.asso[5:6, 1] <- 3
df.asso[c(1, 3, 5), 5] <- "UBI"
df.asso[c(2, 4, 6), 5] <- "TANF"

### All respondents
temp <- prop.test(x = sum(df$asso.UBI), n = sum(!is.na(df$asso.UBI)))
df.asso[1, 2] <- temp$estimate
df.asso[1, 3] <- temp$conf.int[1]
df.asso[1, 4] <- temp$conf.int[2]
temp <- prop.test(x = sum(df$asso.TANF), n = sum(!is.na(df$asso.TANF)))
df.asso[2, 2] <- temp$estimate
df.asso[2, 3] <- temp$conf.int[1]
df.asso[2, 4] <- temp$conf.int[2]

### Conservatives only
temp <- prop.test(x = sum(df[which(df$conservative == 1), ]$asso.UBI), 
                  n = sum(!is.na(df[which(df$conservative == 1), ]$asso.UBI)))
df.asso[3, 2] <- temp$estimate
df.asso[3, 3] <- temp$conf.int[1]
df.asso[3, 4] <- temp$conf.int[2]
temp <- prop.test(x = sum(df[which(df$conservative == 1), ]$asso.TANF), 
                  n = sum(!is.na(df[which(df$conservative == 1), ]$asso.TANF)))
df.asso[4, 2] <- temp$estimate
df.asso[4, 3] <- temp$conf.int[1]
df.asso[4, 4] <- temp$conf.int[2]

### Non-conservatives only
temp <- prop.test(x = sum(df[which(df$conservative == 0), ]$asso.UBI), 
                  n = sum(!is.na(df[which(df$conservative == 0), ]$asso.UBI)))
df.asso[5, 2] <- temp$estimate
df.asso[5, 3] <- temp$conf.int[1]
df.asso[5, 4] <- temp$conf.int[2]
temp <- prop.test(x = sum(df[which(df$conservative == 0), ]$asso.TANF), 
                  n = sum(!is.na(df[which(df$conservative == 0), ]$asso.TANF)))
df.asso[6, 2] <- temp$estimate
df.asso[6, 3] <- temp$conf.int[1]
df.asso[6, 4] <- temp$conf.int[2]

### Visualize the results
df.asso$sample <- factor(df.asso$sample,
                         levels = c(1:3),
                         labels = c("All Respondents",
                                    "Conservatives Only",
                                    "Non-Conservatives Only"))
df.asso$estimate <- df.asso$estimate * 100
df.asso$lwr <- df.asso$lwr * 100
df.asso$upr <- df.asso$upr * 100
fig_asso <- 
  ggplot(data = df.asso, aes(x = sample, y = estimate, color = Policy, shape = Policy)) +
  geom_point(position = position_dodge(.5), size = 3) +
  scale_color_manual(values = c("grey0", "grey60")) +
  scale_shape_manual(values = c(20, 18)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), size = .7, width = 0, position = position_dodge(.5)) +
  xlab("") + 
  ylab("Percentage of Respondents Associating\nthe Policy with Democratic Politicians (%)") +
  theme_bw() +
  theme(text = element_text(color = "black", size = 12, family = "Times"),
        legend.justification = c(1, 1),
        legend.position = c(1, 1),
        legend.background = element_blank(),
        legend.box.background = element_rect(color = "black"),
        legend.key.size = unit(1.5, "line"),
        legend.key.height = unit(0, "cm")) +
  coord_cartesian(ylim = c(0, 100))
ggsave(file = "fig_asso.pdf", fig_asso, width = 5.5, height = 4)